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Cosmic rays acceleration at ultrarelativistic subshocks 
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Abstract 

Wc present the model of cosmic rays acceleration at ultrarelativistic subshocks and confront it with the 
observations of gamma-ray bursts (GRBs) and blazars. We investigate cosmic rays acceleration in shocks 
with Lorentz factors (7) in the range 3 - 40. We show that fluctuations of the magnetic field downstream 
of the shock do not play an important role in the acceleration process. Results of numerical simulations for 
shocks with considered Lorentz factors and perpendicular magnetic field inclinations are presented. We fit 
the derived particle energy spectral index (a) dependence on fluctuations of the magnetic field upstream 
and 7 with a function. 
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1. Introduction 

The acceleration mechanism which operate at ultra- 
relativistic shock fronts was discovered by Bednarz & 
Ostrowski (1998). The mechanism is different from the 
diffusive particle acceleration which is assumed to be sup- 
pressed at superluminal shock fronts (Bell 1978; Drury 
1983). The effect of the magnetic field direction is im- 
portant for ultrarelativistic shocks because all of them 
are superluminal. Medvedev & Loeb (1999) have shown 
that the relativistic two-stream instability naturally gen- 
erate strong magnetic fields which are parallel to the shock 
front. Therefore ultrarelativistic shocks have to be super- 
luminal even if one could imagine an external magnetic 
field with the angle between the upstream field and the 
shock normal smaller than ~ I/7. Ultrarelativistic shocks 
without any mean magnetic fields cannot be considered 
as real physical phenomena by the same reason. 

It have appeared a few papers about a particle accel- 
eration at ultrarelativistic shocks without mean magnetic 
fields (Kirk et al. 2000; Achtcrbcrg et al. 2001; Vietri 
2003; Lemoine & Pelletier 2003) or with subluminal shock 
geometry (Ellison & Double 2002 - parallel shocks) . Their 
acceleration is similar to the diffusive shock acceleration 
but it includes anisotropics in the angular distribution up- 
stream of the shock. However, particles in this accelera- 
tion are able to return to the shock from downstream to 
upstream due to large magnetic field fluctuations down- 
stream of the shock or due to subluminal shock geometry 
as in non-relativistic and mildly relativistic regime. Thus, 
they have failed to understand the actual ultrarelativis- 
tic shock acceleration mechanism because the returning is 
due to small fluctuations of the magnetic field upstream 
of the shock (the needed fluctuations decrease when the 
Lorentz factor of the shock increases) and relativistic ef- 
fects providing a small change of the particle trajectory in 
the mean field upstream of the shock to be large as mea- 
sured downstream. To date, the only numerical calcula- 



tions performed by Bednarz & Ostrowski (1998), Bednarz 
(2000) and in this paper consider the problem of particle 
acceleration in ultrarelativistic shocks. 

It is known for a long time that relativistic shocks occur 
in regions of efficient acceleration of leptons. The acceler- 
ation to non-thermal distributions is observed at hot spots 
of extragalactic radio sources, in blazars, GRBs and pul- 
sar wind nebulae. In order to account for the presence 
of these high energy leptons, some authors have tried to 
find the acceleration mechanism. Begelman & Kirk (1990) 
proposed shock-drift acceleration at relativistic shocks to 
operate at hot spots of extragalactic radio sources. In 
the mechanism, particles are accelerated in a single shock 
crossing by drifting parallel (or anti-parallel) to the elec- 
tric field. Afterwards, Hoshino et al. (1992) described a 
process of shock acceleration of positrons to non-thermal 
distributions devoted to account for the synchrotron ra- 
diation observed in the Crab Nebula and hot spots. In 
the model, the gyrating reflected protons downstream of 
the shock dissipate their energy in the form of collectively 
emitted, left-handed magnctosonic waves which are reso- 
nantly absorbed by the positrons immediately behind the 
ion reflection region. The dynamics of the Weibel instabil- 
ity has recently been simulated by several research groups 
using 3D plasma kinetic code. These simulations confirm 
both the generation of the magnetic field and the particle 
acceleration downstream of the shock. The particle ener- 
gising in Silva et al. (2003) simulations (electron-positron 
plasma) is due to pitch angle scattering in the produced 
magnetic field after saturation. The energy stored in the 
magnetic field is transfered back to the plasma particles. 
Simulations of Frederiksen et al. (2003) have yielded the 
energy transfer from protons to leptons similar to Hoshino 
et al. (1992). In their description, the scattered protons 
create a fluctuating electric field which tends to equilibrate 
the energy between protons and electrons. Nishikawa et 
al. (2003) results suggest that electrons and ions arc ac- 
celerated in the direction transverse to the shock normal 
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only. 

All the described mechanisms suffer from small ener- 
gies the particles are able to acquire. The shock-drift 
acceleration allow for the energy increase of about ten 
times. Leptons can receive the energy from protons which 
is about ten times (~ 7 times in Frederiksen et al. 2003 
does not necessarily depend on 7) above the thermal en- 
ergy downstream of the shock. Silva et al. (2003) simula- 
tions have led to the generation of high-energy tails in the 
distribution function, with energies few times above the 
thermal energy downstream of the shock. In our model 
presented below, we apply some of these mechanisms to 
production of seed particles. 

2. The model of cosmic rays acceleration at ul- 
trarelativistic subshocks 

In the model, we assume an extremely relativistic out- 
flow of electrons and positrons and probably protons and 
nuclei which originates from a central engine. The outflow 
consists of few portions of uniform plasma. These portions 
have different Lorentz factors and relative Lorentz factors 
between some of them reach value from several up to tens. 
Temperatures of these portions are similar to each other. 

In the first phase of the acceleration process, leptons 
and baryonic matter are preaccelerating cither in the 
shock generated by the first slab of the outflow hitting 
the material around the central object or in the shock 
created by two colliding slabs. Depending on 7, seed par- 
ticles are producing in the forward and reverse shock or 
in the forward shock only. 

After the collision the temperature of the resulting 
plasma is higher than temperatures of the two slabs be- 
fore the collision. We expect that the external magnetic 
field influences to the dynamics of the Wcibcl instability 
in such a way that the magnetic field produced by the 
instability is oriented in the form of a toroidal magnetic 
field. It has been suggested that hoop stresses associated 
with toroidal fields wrapped around the extragalactic jet 
will pinch the outfiowing plasma and exert in this way a 
coUimating force on the jet. Following that example, we 
expect that the heating plasma will only expand in the 
direction parallel to its velocity. In result, the collision 
yields a very long and hot jet filled with seed particles 
and toroidal magnetic field. 

We distinguish between dense enough and very thin ma- 
terial of electrons, protons and probably nuclei surround- 
ing the site of the central engine. If the material is very 
thin, then protons or nuclei must be present in the out- 
flow to allow for leptons acceleration thanks to the energy 
transfer from baryons (mechanisms: Hoshino et al. 1992, 
Frederiksen et al. 2003). If the preacceleration takes place 
in two colliding slabs of electron-positron plasma, then 
Silva et al. (2003) mechanism is being the only one at 
work. 

In the second phase, a following slab of plasma is catch- 
ing up with the plasma filled with non-thermal relativis- 
tic particles. These particles arc rapidly accelerating 
in the ultrarelativistic shock acceleration mechanism if 



the Lorentz factor of the formed shock is large enough. 
Afterwards, they flow downstream of the shock and con- 
tinually radiate synchrotron radiation. The maximum en- 
ergy the particles are able to reach is determined by the 
size of the shock or by the value of the magnetic field 
fiuctuations upstream of the shock and 7 or by the accel- 
eration time (it is very short, see Bednarz 2000) and the 
plasma rest frame magnetic and radiative energy density. 
However, only a few per cent of seed particles are involved 
in the process by reason of small efficiency of cosmic ray 
reflections (Bednarz & Ostrowski 1999). 

The presented scenario is capable of explaining many 
features of GRBs and some features of blazars as well. 
The occurrence of the two phases of the acceleration is 
one of them. Ghisellini et al. (2002) have derived from 
spectral properties of blazars that a phase of pre-heating 
and a phase of rapid acceleration leading to a non-thermal 
distribution occur by turns in these objects. In GRBs, the 
first phase is observed in the form of the precursor which 
was detected in a few per cent of GRBs only (Koshut et al. 
1995). The rare occurrence of precursors is best explained 
by small energy of seed particles that strongly depends on 
7- 

In the main phase of the acceleration, the shock slows 
down in consequence of collecting slower material from 
upstream. In the phase we distinguish two sub-phases. 
First, both the reverse and the forward shocks are present 
and the speed of the forward shock is constant. Next, 
the reverse shock disappears on the edge of the back slab 
and the forward shock slows down. The long and hot slab 
of plasma filled with seed particles and toroidal magnetic 
field is much longer than the slab of plasma originating 
from the central engine. In consequence, the reverse shock 
disappears in much shorter time than the forward shock 
crosses the plasma and we can neglect the first sub-phase. 

Further on, we consider that the maximum energy of 
particles in GRBs is determined by the value of magnetic 
field fiuctuations upstream of the shock and 7 (we neglect 
the fraction of particles with large spectral indices) and 
consider constant spectrum of magnetic field fluctuations 
upstream of the shock. As will be shown further, the spec- 
tral index of accelerating particles with momenta larger 
than a threshold momentum pt increases and simultane- 
ously, it remains for the others approximately constant 
as long as 7 decreases. Moreover, pt decreases when the 
shock slows down. The rising phase of the GRB pulse 
is attributed to the energising of seed particles but the 
peak in the light curve is shifted in different energy chan- 
nels. The GRB lag (Kocevski & Liang 2003) arises from 
decreasing pt. When first particles reach pt{t = 0) with 
a K, ai (one should choose a slightly larger than ai - the 
limiting spectral index), the following ones reach pt{t = 0) 
with a> ai because pt{t >0) <pt{t = 0). In result, the ex- 
cess of particles at pt moves towards smaller momenta pro- 
ducing the observed lag. Similarly, the detected spectral 
evolution of main pulses (Crider et al. 1997; Ghirlanda et 
al. 2002) arises due to increasing the energy of particles 
with constant cr w cr; in the rising phase and increasing a 
of particles with momenta larger than pt in the decaying 
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phase. 

Many shapes of main pulses resemble a FRED (fast rise, 
exponential decay). The rising phase of pulses produced 
by shocks with larger initial 7 must yield steeper light 
curves than the phase produced by slower shocks. The 
main reason of this is the reduction in the compression of 
the magnetic field which, in turn, causes the increase of 
the acceleration time (equal to rg/c - Bednarz 2000, rg 
- is the particle gyroradius and c - the speed of light). The 
decaying phase is long because particles are still acceler- 
ating but it can be fast if the shock producing the main 
pulse comes to the region without seed particles before 
particles reach pt. The reverse shock could supply with 
seed particles less efficiently than the forward shock can 
and the light curve will exhibit two peaks when a following 
shock ploughs into such plasma. One could also imagine 
a break in the field of seed particles if two shocks generate 
it independently. 

We propose that short variability timescale pulses, 
which overlap on main pulses, arise from interaction be- 
tween the magnetic field generated by the shock and a 
strong external magnetic field or the toroidal field of the 
hot forward plasma. That interaction can produce fast 
changes in direction of the shock velocity (similar to gyra- 
tion) and, by turns, the sub-pulses because the detected 
flux strongly depends on the angle between the shock nor- 
mal and the direction of observation. 

Most of the particle spectral indices in blazars 
(Ghisellini ct al. 2002) and in GRBs are in good agree- 
ment with the value of 2.18 and larger derived for ultra- 
relativistic shocks in this paper. Suggesting smaller value 
of a (Panaitescu & Kumar 2002) is naturally explained by 
a jet model with continuous energy injection and a above 
2.18 (Bjornsson et al. 2002). The continuous energy in- 
jection is in agreement with the presented model. In some 
of the GRBs events, the last slab (or a group) of plasma 
can be the slowest one. First, the distance between the 
forward and the backward plasma is increasing and the 
forward plasma is collecting interstellar matter and slow- 
ing down. Next, the distance is decreasing and finally, 
the backward plasma hits the forward one and forms a 
shock. We expect that, in many cases, 7 of the shock is 
large enough to allow for the acceleration across the main 
mechanism. The forward material is in the form of long, 
thin jet with the Lorentz factor considerably smaller than 
at the GRB stage. The external magnetic fleld is much 
weaker than the field near the central engine. All the rea- 
sons influence that the detecting radiation is below X-ray 
range and is lasting for long time. The colliding jets are 
responsible for breaks in the afterglow light curves. 

There is an open question if relativistic shocks can 
produce seed protons and nuclei. The simulations of 
Nishikawa et al. (2003) suggest that they could. Heavy 
particles in a cold plasma upstream of the shock are see- 
ing a shock generated by leptons as a discontinuity but 
they are not ultrarelativistic there. We suggest that seed 
baryons could be produced in a mechanism similar to 
shock-drift acceleration in these shocks if leptons taking 
away the energy from baryons downstream of the shock 



remain the fraction of total baryons which end up with 
super-thermal energies. In this picture, the ultrarelativis- 
tic shock acceleration mechanism could produce a major 
part of the cosmic rays. 

3. Simulations 

In simulations we follow the procedure used by Bednarz 

6 Ostrowski (1996) with a hybrid approach used by 
Bednarz & Ostrowski (1998). Monoenergetic seed parti- 
cles are injected at the shock with the same initial weight 
factor. Each particle trajectory is followed using numeri- 
cal computations until the particle escapes through the 
free escape boundary placed far downstream from the 
shock or it reaches the energy larger than the assumed 
upper limit. These particles are replaced with the ones 
arising from splitting the remaining high-weight particles 
in a way to preserve their phase space coordinates, but as- 
cribing only a half or a smaller respective part of the orig- 
inal particle weight to each of the resulting particles. All 
computations are performed in the respective upstream 
or downstream plasma rest frame. Each time the particle 
crosses the shock its momentum is Lorentz transformed 
to the respective plasma rest frame and, in the shock nor- 
mal rest frame, the respective contribution is added to the 
given momentum bin in the particle spectrum. We derive 
the particle trajectories in the mean magnetic field and 
a fluctuating component. The perturbed magnetic field 
represents the traditional picture based on the concept 
of magnetic scattering centres. It is simulated by small 
amplitude particle momentum scattering within a cone 
with angular opening Afi less than the particle anisotropy 
~ 1/7. The particle momentum scattering distribution is 
uniform within the cone. The time (measured in rg/c) 
between scatterings is uniformly distributed from 0.5td,u 
to l.btd^u (td and i„ are mean times between scatterings 
downstream and upstream respectively). After crossing 
the shock, the remaining time the particle should follow 
in the mean field before the next scattering is multiplied 
or divided by td/tu- Particle energy is preserved at each 
scattering in the plasma rest frame. For the considered 
continuous injection after initial time, the energy cut-off 
of the formed spectrum shifts toward higher energies with 
time. The resulting spectra allow one to fit spectral in- 
dices and derive acceleration time. The downstream mag- 
netic field, its inclination to the shock normal and shock 
speed arc derived for the relativistic shock with the com- 
pression Rh = "iUi {Ui is the shock speed) for the hot 
plasma. In the simulations presented in Fig. 2. we applied 
compression Rc obtained with the formula of Heavens & 
Drury (1988) for a cold proton-electron plasma. 

In order to estimate the correctness of Rh for slow 
shocks, we followed the computations of Heavens & Drury 
(1988). The bottom panel of Fig. 1 shows the difference 
between R (the real compression) and Rh as a function of 

7 for some cases of the upstream fluid (Tp, Tg - the up- 
stream temperature of proton-electron, electron-positron 
plasma in units of 5.4 ■ lO^^if , 5.9 • lO^i^ respectively). 
The top panel shows the relation between 7 and the down- 
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Fig. 1. The downstream temperature and the difference be- 
tween the compression and the compression of the hot plasma 
limit as a function of the shock Lorentz factor. The curves cor- 
respond to different upstream temperatures and are marked 
with Tp (solid curves) or Tg (dotted curves). In the upper 
panel, Tp = Te = 10'^° is close to Tp = 0.01. 

stream temperature (T^ - in the same unit as appropriate 
Tp_e). We can draw a conclusion that 7 of the shock hitting 
a cold plasma must be larger than 6 during the preacceler- 
ation so that Rh might be adequate for shocks with 7 < 5 
in the main phase. 

4. Some features of the acceleration mechanism 

We have performed simulations in the absence of mag- 
netic field fluctuations upstream and for a few cases of 
the fluctuations downstream measured by a parameter 
A = log{^), where k± and K|| are perpendicular and par- 
allel diffusion coefficients respectively. We have chosen a 
few values of 7 and the upstream magnetic field inclina- 
tion with respect to the shock normal = 10°, 30°, 90°. 
We exhibit the outcome of the simulations in Fig. 2 . 
The large value of ct for -0 = 10° and 7 = 8 arises due 
to escaping a large fraction of particles far upstream from 
the shock. The a dependence on 1/) or 7 is present for 
large value of a only. We have obtained a = 2.5; 3.5; 5 for 
A = 0; —0.16; —0.67 respectively. 

As was shown above, medium- value magnetic field fluc- 
tuations downstream can be neglected in the accelera- 
tion. We suspect that they are really small if we take 



into account very high linear polarization reported in 
RHESSI observations of GRB021206 (Coburn & Boggs 
2003). Moreover, the particles do not interact resonantly 
with magnetic field perturbations upstream since the 
time between shock crossings downstream-upstream and 
upstream-downstream is a small part of rg/c. Therefore, 
our scattering model is the best because it reproduces the 
non-resonant interactions excellently and independently 
of magnetic field fluctuations spectra. The model follows 
the relation ^ = -r-, hi — vs, t is the mean time when 

K|| l+const{ujT)^ ' 

the particle momentum (p) direction shifts at an angle 
smaller than 20° because of scatterings and uj is the gyra- 
tion frequency. Let us introduce the ordinary assumption 
that each fluctuation shifts the direction at a small angle 
^ 1/p. As a result, t ^ and subsequently, ^ = i^Qpi 
{Q results from fluctuations spectra) . For large p wc get 
X = —2logio{Qp)- Eq. (1) yields that A upstream increases 
{p decreases) with decreasing 7 and constant a. That is 
the reason for the lag in GRBs light curves. 

The angular distribution of particles at the shock in 
the downstream plasma rest frame presented by Bednarz 
& Ostrowski (1998) peaks at fi = —U2, U2 is the shock 
speed in the downstream rest frame and fi is the cosines 
of the angle between the shock normal and the particle 
momentum. The distribution presented by other authors 
peaks at /i ~ but they did not follow full 3D simula- 
tions (or calculations) with the mean field and we cannot 
compare the results because it is difficult for us to mimic 
the reduced situations. We expect the peak at —U2 arises 
from keeping a fraction of particles close to the shock. 
The particles contribute to the distribution considerably 
because their velocities perpendicular to the shock plane 
approach the shock speed. They can oscillate between up- 
stream and downstream especially, if scatterings upstream 
are close to the shock plane. In Fig. 3, we present the sim- 
ulations for 7 = 27, tp = 90°, different td^u (the diffusion 
is the same because Ail'^ /td^e = const) and A = 0. The 
distribution function narrows with decreasing td.u but a 
does not change. 

The explanation of the shift in the distribution function 
presented by Ostrowski & Bednarz (2002) is false^. It 
is obvious that the absolute value of the inverse of the 
relative particle velocity with respect to the shock is the 
particle weight. One obtains unlike distribution if one 
applies the velocity in the plasma rest frame. 

We obtained the limiting spectral index for A = up- 
stream and downstream, ip = 90° and 7 from 60 to 120 
(33 points). It is equal to crim{j > 60) = 2.23 ±0.03. 

5. The relation between cr, 7 and A 

In the simulations below, we have applied td = 0.1 and 
A57 = 10~^ downstream (e^ = 0, homogeneous magnetic 
field) and tp = 90°. First, we have found the limiting spec- 
tral index for large 7 (7 from 60 to 100 - 58 points. A = 
upstream): cr/o(7 > 60) = 2.18 ±0.03. Next, for the same 

^ That is not mine. M. Ostrowski did not consult with me the 
context of the paper. 



No. 



14- 
12 ■ 
10- 
8- 
6- 
4- 

10 

8- 
6 
4 
2 - 
10 
8- 
6- 



A, = -0.67 



H h 

X = -0.16 



G-o 10 

30 

"■■■« 90 



H h 

X = 0.00 



27 
Y 



81 



243 



Fig. 2. Simulated spectral indices as a function of the shock 
Lorentz factor. Three values of the magnetic field inclina- 
tion upstream of the shock are applied (i/i = 10°, 30°, 90°). 
Magnetic field fluctuations upstream are not present. Results 
for given downstream magnetic field perturbations A are pre- 
sented separately. 
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Fig. 3. Simulated particle angular distributions at the shock, 
as measured in the downstream plasma rest frame, for 7 = 27, 
4' = 90° and A = (downstream and upstream). The distribu- 
tion function narrows when „ decreases but a is constant. 




Fig. 4. The difference between the limiting spectral index 
obtained from simulations (cris) and cr;. 

conditions {ip = 90° and ~ downstream), our simu- 
lations have covered two 7 ranges: 1) from 2 to 7 every 
0.1; 2) from 2 to 40 every 1 and A upstream from to 
-4.3. Then, we have been able to obtain cr; for 7 > 3: 
(Ji = 5.18 — Rh (see Fig. 4). We have continued analysis of 
the data. Let us introduce two parameters Tm — Rh/^ and 
rd = Rh — 2. We have fitted a/ai — r^;^ as a function of 
Ap = 3.36 - 1.94^0310(7) - 0.3r<j 
and 3 < 7 < 40: 



X/rd for o- < 7, Ap > 0.55 



13.1Ap)(l- 



0.147^ 



= 4.6 • 10-3 ± 



2-10' 



(1) 



The fitting is correct if Ap > but the error is too large and 
there are systematic effects if cr ^ cr/ . We have replaced (1) 
with the linear function from cr « 1.03cre [u^ is the limiting 
cr derived from (1)) to cr/ to improve eq. (1). We have fixed 
the intersection of the two curves BiXi A-n 



A 



2.86 - 0.2657+ 0.01572 - 4.02 
MO^V- 



Ap(cr w 1.03cr, 

7^+ 



(2) 



{{^)c-l)/K, (3) 
r3„2)/(2.3Ap,) de- 



The linear function is 

cr = cr/(aA+ 1) ± Acr, 

where the select error Acr = cr/((^) 
pends on 7 only at the intersection (^)c — 4.6 • + 
13.1Ap,)(l - 23MZ) + r3„2, K = TdiiM - 1.94^0510(7) - 
Q.ird — Ape) and below. An example of a curve derived 
from (1), (2) and (3) for 7 = 32 is presented in Fig. 5 and 
a few curves for different 7 are shown in Fig. 6. The rela- 
tion between cr, 7 and A is necessary to model GRBs light 
curves. 
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The dotted lines represent the fit. 
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